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' Abstract 

bJO' We consider a spherically symmetric, static system of a Dirac particle interacting 

^ i with classical gravity and an SU{2) Yang-Mills field. The corresponding Einstein- 

' Dirac- Yang/Mills equations are derived. Using numerical methods, we find different 

types of soliton-like solutions of these equations and discuss their properties. Some of 
these solutions are stable even for arbitrarily weak gravitational coupling. 
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! 1 Introduction 

The coupling of gravity to other classical force fields and to quantum mechanical particles 
\ has led to many interesting solutions of Einstein's equations and has given some insight 

into the nature of the nonlinear interactions. The first such examples are the Bartnik- 
McKinnon (BM) solutions of the Einstein- Yang/Mills (EYM) equations Q. For these 
qh| solutions, the repulsive Yang-Mills force compensates the attractive gravitational force; 

unfortunately, these solutions are unstable Q. If, on the other hand, one considers quan- 
bi)' tum mechanical Dirac particles, the gravitational attraction is balanced by the repulsion 

due to the Heisenberg Uncertainty Principle, and this leads to stable bound states of the 
resulting Einstein-Dirac (ED) system Q. However, a pure ED system is somewhat ar- 
tificial, because physical Dirac particles also exhibit electroweak and strong interactions, 
5^ \ which in all realistic situations are much stronger than gravity. Thus the question arises 

if Dirac particles in a gravitational field still form bound states if an additional strong 
coupling to a non-abelian Yang-Mills (YM) field is taken into account (the case of an 
abelian gauge field was considered in [Q]). Related questions are, do the BM solutions be- 
come stable if one includes Dirac particles, and which physical effects does the nonlinear 
coupling in the Einstein-Dirac- Yang/Mills (EDYM) equations lead to? 

In order to address these questions, we consider here a static, spherically symmetric 
EDYM system of one Dirac particle in a gravitational field and an SU (2) Yang-Mills field. 
In this system, the spinors couple only to the magnetic component of the YM field, and we 
thus obtain a consistent ansatz by setting the electric component identically equal to zero. 
In the limit of weak coupling of the spinors, our system goes over to the EYM system Q. 
In contrast to the two-particle singlet state studied in , we consider here only one Dirac 
particle (this becomes possible because the inclusion of the SU{2) Yang-Mills field changes 
the representation of the rotation group on the spinors; see Section 0). Thus one cannot 
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recover exactly the ED system Q , but the hmit of a weak Yang-Mihs field yields equations 
which are closely related to the ED equations of the two-particle singlet. 

By numerically seeking bound states of the EDYM system, we find a surprisingly rich 
solution structure. First of all, we find solutions where the Dirac particle is bound by the 
gravitational attraction, and where the Dirac particle also generates a YM field. Stable 
solutions of this type exist also in the physically realistic situation of weak gravitational 
and strong YM coupling. This result shows that the magnetic component of the YM 
field, which usually has a repulsive effect (like e.g. for the BM solutions) cannot prevent 
Dirac particles from forming stable bound states. We also find other types of solutions 
where the binding comes about through the nonlinear interaction of all three fields. These 
solutions have stable and unstable branches, whereby the stable solutions are found for 
weak gravitational coupling, provided that the YM coupling is sufficiently strong, but not 
too strong. Finally, we study the relation between these solutions and the BM solutions. 
We find one-parameter families of solutions which join the BM ground state with our new 
solutions. This shows that the BM ground state can be made stable by the inclusion of 
a Dirac particle, but only if the coupling to the spinors is sufficiently strong. The first 
excited BM state, on the other hand, cannot be joined with our new solutions. Namely, 
perturbing this state by an additional Dirac particle yields a separate unstable branch of 
EDYM solutions. 

The plan of the paper is as follows. In Section ^ we derive the SU (2)-EDYM equations. 
In Section ^ we obtain a limiting system constructed by letting the gravitational constant 
tend to zero and letting the YM coupling constant tend to infinity. We find numerical 
solutions of this system and discuss their properties. In the last section, we consider 
solutions of the full EDYM equations, obtained by tracing the solutions of our limiting 
system and the BM solutions while continuously varying the coupling constants. 



2 Derivation of the Equations 

The general EDYM equations are obtained by variation over Lorentzian metrics gij, YM 
connections A, and Dirac wave functions of the action 

S = [ f-^— R + ^iG- m)^ - Tr{Fij r^)) V^d^ d^x , (2.1) 

where R is scalar curvature, G is the Dirac operator (which depends on the gravitational 
and YM fields), Fij is the YM field tensor, and where the trace is taken over the YM 
index. The gravitational and YM coupling constants are denoted by k and e, respectively. 



The appearance of the factor e"^ in {2A.) requires a brief explanation. In contrast to the 
usual form of the gauge-covariant derivative Dj = dj — ieAj, we use here the convention 
Dj = dj — iAj (this makes it possible to work with the particularly convenient form of 
the gauge potentials used in [l|). Our convention is obtained from the usual one by the 
transformation Aj e^^ Aj. Under this transformation, the field strength tensor behaves 



like Fij — > Fij, and this gives rise to the factor e"^ in (2.1). 

In this paper, we shall study a particular EDYM system, which is obtained as follows. 
First of all, we consider a spherically symmetric, static metric in polar coordinates, 

ds^ = -^dt" - -^dr' - rU^^ - sin^^? ^ ^ (2.2) 
T[rY Ayr) 
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with positive functions A and T. The Einstein tensor corresponding to this metric is given 
in 1^. Moreover, as in we restrict attention to the magnetic component of an SU{2) 
Yang-Mills field and choose for the YM potential the ansatz 

A = w{r) d-d + {cos-dr^ + w{r) sin i^T^)dip (2.3) 

with a real function w, where r = ^ <? is the standard basis of su{2) {a are the Pauli 
matrices). The YM current j and energy-momentum tensor Tj = Tr {F^"'Fja — 
1 F°-''Fab 5j ) are computed to be 

' ^ J-TTT^^ - K 11^ + ^ csc(7?) 



47re2 \^ 2r2 4r2 T 2r'^ ) \ ^ ' dip 



47re2 V '''^ ''"^ 



= Ti = , (2.6) 

and all other components vanish. When coupled to the Dirac spinors, the YM potential 
( ^.31 ) has the disadvantage that it depends on and ip, in a way which makes it difficult 
to separate variables in the Dirac equation. To remedy this, we perform the SU (2) gauge 
transformation Aj UAjU~^ + iU{djU~^) with 



U{'d,(p) = ex.p (^—iip exp ^— i(?? -|- vr) r'^^ exp ^- 

The resulting YM potential is 

A = {w - 1) r sin 1^ {t'P - t"^ dip) , (2.7) 

where we use the following "polar" linear combinations of the r matrices, 

= COS"!? -|- siwd cosp + siwd sin 95 
t"^ = - (^—T^ siwd -|- cos-d cos(/? -|- cos-!? sin 99 



t'^ = (— sin 99 -1- cos 93 

rsmw V 



By minimally coupling the SU{2) potential ( p.TD to the Dirac operator in the gravitational 
field Eq. (2.23)], we obtain the Dirac operator 

G = iT ^'dt + Y (iVAdr + -{Va-i)--!-Va^] + i^d^ + i^d^ 

V r 2 1/ 

+ -{w-l){jT-Yr'^)r'- , (2.9) 
r 



where ('y^)j=t^r,&,ip are, in analogy to (|2.8D , the Dirac matrices of Minkowski space in polar 
coordinates, where we work in the Dirac representation, 
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Notice that the Dirac operator (|2.9D acts on 8-component wave functions; this is because 
the additional YM index doubles the number of components compared to usual Dirac 
spinors. More precisely, it is convenient to regard the wave functions as sections of 

= ^up/down ^Lge/small ® ^YM ) (^-H) 

where <Cup/down describes the two spin orientations, ^Cfargc/smaii corresponds to the upper 
and lower components of the Dirac spinor (i.e., usual Dirac spinors are sections of = 
<Cun/down ® ^Lee/small)' ^ud (C^^ is acted upon by the SU{2) gauge group. For clarity, 



we shall refer to the factors in ( 2.11 ) by separate indices, i.e. we write a wave function ^' 
as (^°"'')a,n,a=i,2, where a, n, and a label the components of (^1^/^^^^, (CLge/smaii' ^"^^ 
(DyM) respectively. Thus the operators r act on the index a, the spin operators S are given 
by 5 = ^ (T acting on the Greek indices, and 7* coincides with the matrix 7* = diag(l, —1) 
acting on the index u, i.e. 



7* ^''^ 



It is apparent in ( |2.9D that the Dirac operator commutes with the three operators 

J = L + S + T , (2.12) 

where L is angular momentum. It is very convenient to regard the operators J as the total 
angular momentum operators of the system. Since the total angular momentum operators 
are the infinitesimal generators of rotations (as explained for angular momentum in |5|, par. 



26]), we can then interpret ( 2.12 ) as saying that the inclusion of the YM field influences 
the representation of the rotation group on the spinors. The Dirac operator is invariant 
under this group representation, because the operators J commute with G; this makes 
spherical symmetry of the Dirac operator manifest. 

Since ( 2.12| ) coincides with the formula for the addition of angular momentum L to 



two spin-^-operators S and r, it is clear that the operators J have integer eigenvalues. 
Thus we can expect that the operator has a nontrivial kernel. In this case, the simplest 
spherically symmetric ansatz for the Dirac particles would be to take one Dirac particle 
whose wave function is in the kernel of J^. We now work out this ansatz in detail, whereby 
we consider J as operators on the spinors <I>°°(t?, 93) on S'^ (i.e. the ^'^'^ are sections of 
(C^ = tCup/down ^ym)- Adding the two spin operators S and r, we can decompose 
^up/down ® ^YM ^he direct sum of one state of total spin zero and three states of 
total spin one (see [||, par. 31]; these states are usually called the singlet and triplet states, 
respectively) . By subsequently adding the angular momentum L according to the standard 
rules for the addition of angular momentum ||5|, par. 31], one sees that the operator has 
indeed a nontrivial kernel. More precisely, the kernel of is two-dimensional, spanned 
by two vectors $0 and $1 with angular momentum zero and one, respectively. The state 
'I'o is (up to a phase) uniquely characterized by the conditions 

L^>o = = (5^ + f)$o and H^oll^a = 1 . (2.13) 

Using (2.13), we can write <I>i as 

$1 = 25"^ $0 = -2t''^o • (2.14) 
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Namely, representing S'^ and in the form 

= xs and = xt , 

and using the standard commutation relations between the components of L, x, and S 
(see ||5|, pars. 26 and 54]), we obtain that 

JcDi = 2 [J, 5n $0 = 2 [L, cDq + 2 [5^ + f, cDq 

= -2ix a5^>o + 2ix a5$o = (2.15) 

(where A is the wedge product in IR'^), and 

\\<^ifs2 = I <2S''^o,2S''^o>duj = ||$o|||2 = 1 . 

One can verify directly that has angular momentum one; namely 

= 2L^ S'' = 2L[L, S'']^Q = -2iL{x hS)^Q 
= -2i [L, (x A 5)] $0 = 4 {xS) $0 = l{l + 1) $1 

with I = 1. Furthermore, using the fact that (25"^)^ = 1 = (2r'')^ and 5"^ = | = r^, we 
obtain that 

$0 = 25^ $1 = -2r''$i (2.16) 

St^q ^'^=^ -SS^o = -|$o (2.17) 

St^^ = ^{{S + rf - S2 - T^) $1 S i (l2 - - t') <I>i = ^ 1>i (2.18) 
SL$o = (2.19) 
SL<^i = ^((5 + L)2-52-L2)$i = i (r2 - 52 - l2) $1 = -^.^ . (2.20) 

Finally, it is useful to observe that (cf. [^, equation (3.3)]) 

S^d^ + S'^d^ = --S^iSL) . (2.21) 
r 

Using the relations (^l3|)-(g^, one can easily compute the Dirac operator (|2.9| ) on the 
invariant subspace = 0. It turns out that we obtain a consistent ansatz for the Dirac 
wave function by setting 



^-"'^(t, r, 1?, 99) = e-'"* ^^-^ (a(r) $^^^(1?, <^) (^,,1 + i/3(r) $r (^^, <2) (2.22) 

with real functions a and /3, where u; > is the energy of the Dirac particle, and is the 
Kronecker delta. For this ansatz, the Dirac equation reduces to the system of ODEs 

^/Aa' = -a-{m + ojT)f3 (2.23) 



r 



= {-m + ujT)a (3 . (2.24) 

r 
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The Dirac current j and Dirac energy-momentum tensor Tj/. = Re (^'G(j-Dfc)^) corre- 
sponding to the ansatz ( |2.22| ) are obtained by a straightforward computation similar to 
that in [01. The result is 



2T 



a 9 



d 



Oif^\'^ ^ + csc(t?) — 



dip 



Tl 



T 

— w aP 



and all other components vanish. The normalization condition for the spinors is (as in 

T 



JO 



(2.25) 



By substituting the formulas for the YM current and energy-momentum tensor into 
the Einstein and YM equations, we obtain the following system of ODEs, 



rA' 

r 



12 



2rA 



r^A w" 



.l + A + + 2Kmr(a2_/?2) 

-4k — w ap — — r- A w 



-{l-w'^)w + e'^rTaP 



oA'T-2AT' , 

r w 

2T 



(2.26) 

) 

(2.27) 
(2.28) 



The Einstein equations are ( 2.26 ) and ( |2.27| ), whereas ( 2.28 ) is the YM equation. Our 
EDYM system is given by the five ODEs (|2?^ ), ([2?24|) , ( |2?26| )-( |2^ , together with the 
normalization condition ( p. 25 ). 

We are interested here in bound states of the Dirac particles. Thus we want to find 
particle-like solutions of our EDYM system, i.e. solutions which are smooth and tend to 
the vacuum solution as r ^ oo. According to the explicit formulas (2.4)-(2.6), the energy- 
momentum tensor of the YM field is regular at r = only when 11(^(0)1 = 1 and w'{Q) = 0. 
Using the remaining gauge freedom, we can assume that w{Q) = 1, and thus 



w{r) 



1 - - r^ + 0{r^) 



(2.29) 



with a real parameter A. Using this result, a local Taylor expansion of the Einstein and 
Dirac equations around r = yields (just as in S) that 



a(r) 
A{r) 



air + 0{r'^) 
1 + 0(r2) 



^ (wro-m)air2 + 0{r^) 



T{r) = To + 0{r^) 



(2.30) 
(2.31) 



with two parameters ai and Tq > 0. Using linearity of the Dirac equation, we can always 
assume that ai > 0. Furthermore, we demand that our solution has finite ADM mass. 



p := lim — (1—A(r)) < oo 

r^oo 2k 



(2.32) 
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and goes asymptotically to the vacuum solution, 

lim r(r) = 1 , lim (w(r),w'(r)) = (±1,0) , lim (a(r), /3(r)) = (0,0) 

r — *-oo r — >oo r — *-oo 

3 The Reciprocal Coupling Limit 



(2.33) 



Under all realistic conditions, the coupling of Dirac particles to the YM field (describing 
the weak or strong interactions) is much stronger than the coupling to the gravitational 
field. Thus we are particularly interested in the case of weak gravitational coupling. In 
preparation, it is instructive to briefly consider the case without gravitation. In this limit, 
the Dirac equations read 

w 

a' = -a - {m + uj)f3 (3.1) 
r 

w 

(3' = (-m + uj) a (3 . (3.2) 

r 

For large r, these equations go over to a linear system of ODEs with constant coefficients, 
and the sign of m—u determines whether the solutions of these equations behave oscillatory 
or exponentially. The normalization condition ( p. 25 ) excludes the oscillatory case (as in |^ 



Section 5]) and thus m — a; > . In the case m — a; = 0, the /3-equation is independent of 
a, and the boundary conditions ( 2.30| ) imply that /3 = 0. As a consequence, (g) reduces 



to the homogeneous YM equation 

r'^w" = -{l-w'^)w . (3.3) 

It is well-known [0] that the only solution to this equation satisfying the boundary con- 
ditions ( |2.29| ),( p.33| ) is the trivial solution w = 1. But then the a-equation simplifies 
to 

, 1 

a = — a , 
r 



whose solution a = air violates the normalization condition ( 2.25| ). In the case m — uj > 0, 



on the other hand, the local Taylor expansion (2.30) yields that the (a, /3)-curve lies 



for small r in the fourth quadrant, i.e. /3(r) < < a{r) for small r. Using the Dirac 
equations ( |3.1D ,( p^ ), one sees that the fourth quadrant is an invariant region, and thus 
(5{r) < < a{r) for all r. But in the fourth quadrant, both a{r) and — /3(r) are increasing 
for large r (as one sees in ( |3.lD ,( ^^) ta king into account that u;/r^Oforr^oo), and 
thus the normalization condition (|2.25D will again be violated. 

These considerations show that the gravitational field is essential for the formation of 
bound states. Nevertheless, for arbitrarily weak gravitational coupling, we can hope to 
find bound states. It is even conceivable that these bound state solutions might have a 
well-defined limit when the gravitational coupling tends to zero, if we let the YM coupling 
go to infinity at the same time. Our idea is that this limiting case might yield a system of 
equations which is simpler than the full EDYM system, and can thus serve as a physically 
interesting starting point for the analysis of the coupled interactions described by the 
EDYM equations. Expressed in dimensionless quantities, we shall thus consider the limits 

w?Hi — > and — > oo . (3.4) 
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Let us determine how the quantities of our EDYM system should behave in this hmit. 
Since we are considering weak gravitational coupling, it is clear that the metric will be 
close to the Minkowski metric, i.e. A ~ 1 and T ~ 1. Furthermore, the YM potential 
w should have a finite limit. Similar to our flat space consideration at the beginning of 
this section, one sees that the normalization condition ( |2.25| ) can be satisfied only if the 
function m — ujT{r) changes sign, and thus to ^ m (but both m and cv may go to zero or 
infinity in the limit ( p.4D). Putting this information together, we conclude that the Dirac 
equations ( 2.23| ) and ( |2.24|) have a meaningful limit only under the assumptions that a 
converges and that 



m (3{r) /3(r) 



(T(r) 



m {b 



m) 



E 



(3.5) 



with two real functions (5, (p and a real parameter E. Multiplying (|2.2^ ) with m and 
taking the limits (p.5|) as well as A, T — > 1, the Dirac equations become 



w 



a 



a 



2/3 



{E + Lp)a 



w 



(3.6) 
(3.7) 



We next consider the YM equation (|2.28D . The last term in ( ^.28 ) drops out in the limit 
of weak gravitational coupling ( |3.4[ ). The second summand converges only under the 
assumption that 

> q (3.8) 

m 

with q a real parameter, playing the role of an "effective" coupling constant. Together 
with (^.41), this implies that m ^ cc. The YM equations thus have the limit 

,.2\ 



2 // 

r w 



-{l-w^)w + qrap . (3.9) 

In order to get a well-defined and non-trivial limit of the Einstein equations ( |2.26| ),( p.27D , 
we need to assume that the parameter m^n has a finite, non-zero limit. Since this param- 
eter has the dimension of inverse length, we can arrange by a scaling of our coordinates 
that 

m^K 1 . (3.10) 
We differentiate the T-equation ( 2.27 ) with respect to r and substitute ( 2.26| ). Taking the 



limits (3^) and ( 3.10 ), a straightforward calculation yields the equation 

„2 



-a 



(3.11) 



where A 



~'^dr{r^dr) is the radial Laplacian in Euclidean IR^. Indeed, this equation 
can be regarded as Newton's equation with the Newtonian potential (p. Thus our limiting 
case ( p.ll ) for the gravitational field corresponds to taking the Newtonian limit. Finally, 
the normalization condition (2.25) reduces to 



a(r)^ dr 



The boundary conditions ( 2.29 )-( 2.33D are transformed into 



w{r) 

a{r) 
ip{r) 



air + ©(r^) 



lim wir) = ±1 

/5(r) = 0{r^) 
lim ip{r) < oo 



(3.12) 

(3.13) 

(3.14) 
(3.15) 
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with the three parameters A, ai, and (fQ. We point out that the Umiting system contains 
only one couphng constant q. According to ( |3.§| ) and ( 3.10| ), q is in dimensionless form 
given by 

q . (3.16) 



2 2 

emu 



Hence in dimensionless quantities, our limit ( p.4| ) describes the situation where the grav- 
itational coupling goes to zero, while the YM coupling constant goes to infinity like 
~ {m'^K)~^. Therefore, we call our limiting case the reciprocal coupling limit. The 
reciprocal coupling system is given by the equations ( p.6D , ( |3.7| ) , ( |3.9| ) , a nd (|3.11[) to gether 
with the normalization condition ( 3.12 ) and the boundary conditions ( p.l3D -( p.l5D . Ac- 
cording to (p.5|), the parameter E coincides up to a scaling factor with uj — ni, and thus 
has the interpretation as the (properly scaled) energy of the Dirac particle. As in Newto- 
nian mechanics, the potential (p is determined only up to a constant fj, £ IR; namely, the 
reciprocal limit equations are invariant under the transformation 



ip + fi 



E 



E 



1^ 



(3.17) 



Let us consider how the ADM mass behaves in the reciprocal coupling limit. First of 
all, we can write the quotient p/m as 



— = lim 

m r^oo 2Km 



[1-A{r)) = - r (l - Air)))' 
m Jo \2k ) 



dr 



After substituting the j4-equation ( |2.26| ) , we can take the limits 
that 



P_ 

m 



a(r)^ dr 



(3.12) 



1 



and (|3.5|) and obtain 



(3.18) 



Thus the ADM mass coincides with the rest mass of the Dirac particle; this shows that 
the total binding energy B := p — m goes to zero in our limit. Indeed, the behavior of the 
total binding energy can be described in more detail as follows. For a solution of the full 
EDYM system, we can write the binding energy using the normalization condition (|2.25| ) 
as 



B 



-(1 

2k ^ 



A) 



T 

m [a^ + B^) —= ] dr 



We again substitute the A-equation ( ^.26|) and obtain 

,2^2 I 



B 




Aw''^ + {wT^/A - m) — = (a^ + /32) j dr 



(3.19) 



According to ( p. 16 ), it is obvious that the first two summands in ( p. 191 ) have a finite 
limit after dividing by rn^K. In order to treat the last summand, we first multiply the 
T-equation (^37|) with and take the limits (U), (jsj), (|3.10D , 



{A - 1) 



2r (f' 



Using again ( 3.10D , (|3.5D , and u; — > m, T — > 1, we obtain that 



{ujT^/~A — m) 



muo 



{^/A-l)T + 



m(u)T 



m) 



rip' + {ip + E) 
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From this we conclude that the binding energy 
hmit; more precisely 

S := 4- - r (^^-^—^ + + ^'iE + ^ + r^')]dr . (3.20) 

m^K Jq \2q r-^ q J 



3.19) divided by m k has a meaningful 



We now describe our method for constructing numerical solutions of our reciprocal 
limit system. Since it is difficult to take into account the integral condition ( 3.12| ) in the 



numerics, we discard this condition for the construction of the solution; it will be taken 
care of later via a rescaling technique (see ( 3.341 ), ( |3.35| )). This rescaling method requires 



only that the normalization integral be finite, 

POO 

< := / a{r f dr < oo . (3.21) 



According to (|3.6|) , (3.7) and ( ^.13]) , (|3.15|) , the behavior of the Dirac spinors at infinity is 



either oscillatory or exponential. As a consequence, the normalization integral in ( p. 21 ) will 



be finite only if a(r) tends to zero for r — > oo. Furthermore, we can use the transformation 
( |3.17|) to set ip{0) = 0. Hence in the first construction step, we want to find solutions of 
the modified system 

w 

a' = -a - 2(3 (3.22) 
r 

w " 

p' = {E + ip)a (3 (3.23) 

r 

r^w" = -{l-w^)w + qraP (3.24) 
r"^ A^{r) = -a^ (3.25) 

with the following conditions at the origin, 

w{r) = l-Xr'^ + 0{r^) , a{r) = air + Oir"^) (3.26) 

y,(r) = 0(r2) , /?(r) = 0{r^) , (3.27) 

together with the conditions at infinity 

lim w(r) = ±1 , lim a(r) = (3.28) 

r^oo r^oo 

\(Poo\ '■= I lim ^(r)\ < OO . (3.29) 

r — *oo 

For any given value of the coupling constant q, we thus have two free parameters A and ai 
to characterize the solutions near the origin r = 0. Each solution has a unique extension 
to larger values of r. Asymptotically for r ^ oo, we must satisfy the two conditions ( |3.28| ). 
Thus we have as many free parameters as asymptotic conditions, and we therefore expect 
for fixed q a discrete set of solutions satisfying ( p.26|) , ( |3.27| ), and ( 3.2g| ). For each solution, 



we must then verify that the conditions ( p.21[ ) and ( |3.29D are also satisfied. 

For the construction of numerical solutions, we enhanced the technique used in |^ ^ to 
a two-parameter shooting method. Since two-parameter problems are considerably more 
difficult than one-parameter problems, we describe the method in some detail. For clarity, 
we first consider the simplified situation where a{r) and w{r) have prescribed boundary 
values for a given finite r = ri. In this case, one can apply the standard multi-parameter 



10 



shooting method as e.g. described in P|. More precisely, one can for given initial data 
compute a(ri) and {w{ri),w'{ri)) numerically, compare with the prescribed boundary 
conditions, and iteratively adjust the two free parameters A and ai until the boundary 
conditions are satisfied to sufficient accuracy. In our case, the situation is more difficult 
because we have boundary values not for finite r = ri, but for r = oo. In order to deal 
with this problem, we first choose a finite ri. Using an ansatz for the asymptotic form of 
the solution {a,P,w,if) at infinity, we approximately compute a{ri) and {w{ri),w'{ri)) 
and derive conditions between these functions. Taking these conditions as the boundary 
conditions at r = ri , we can apply the two-parameter shooting method on the finite interval 
(0, ri] as described above. The so-obtained solution on (0, ri] gives, in combination with 
the asymptotic formulas on (ri,oo), an approximate solution for all r > 0. Since our 
asymptotic description becomes precise only in the limit ri oo, we must, in order to 
attain the desired accuracy, choose ri sufficiently large. In order to ensure that ri is 
appropriately increased during the computation, we modified the two-parameter shooting 
method in such a way that both the initial data and ri are adjusted in each iteration 
step. The iteration is stopped when the numerics has stabilized and the accuracy no 
longer improves. This modified shooting method was implemented in the Mathematica 
programming language using the standard ODE solver with a working precision of 16 
digits. The initial data is adjusted in the iteration with a secant method, and the step 
size for incrementing ri is determined from the relative error of the numerical solution at 
the upper boundary ri. After the iteration has been stopped and a numerical solution 
has been found, our program slightly changes the initial data and searches for a nearby 
solution. In this way, we can automatically trace a one-parameter family of solutions. 
Finally, we explain our method for describing the asymptotic behavior of the solutions 
at infinity. According to the asymptotics of the solutions of the ED and EYM equations 
[l|, we can expect that the spinors a and (3 will decay exponentially fast at infinity, 
whereas the potentials and w^r) for r — > cxd will behave like rational functions. 

Therefore it is a reasonable asymptotic approximation to set a and (3 to zero for r > ri. 
In this approximation, the potential ip is harmonic according to (p. 251 ). The YM equation 
( |3.24| ), on the other hand, reduces to the vacuum YM equation ( |3.3D . In the new variable 
u = logr, this equation becomes autonomous; namely Q 

dlw-duw = -{l-w'^)w . (3.30) 

This autonomous equation allows us to derive boundary conditions for w as follows. We 
set 

X = w and y = duW . (3.31) 
Then the YM orbits in the (x, y) plane are described by the following differential equation, 

y'ix) = ^'^'^ 1 - ^i^^ . (3.32) 

duW y 



According to the boundary conditions ( p. 28 ) and the differential equation ( 3.30| ), the 



variables x and y must behave in the limit r — oo like either a;^l,y\Oorx^— 1, 
y y 0. In both of these cases, there is a unique YM orbit y{x), which can be easily 
calculated numerically by integrating ( 3.32] ). By transforming ( |3.31D back to the variable 



r, we obtain the following mixed boundary conditions for w{r) at r = ri, 

w'{ri) = ^y{w{ri)) . (3.33) 
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We next describe our rescaling method needed to arrange the normahzation condition 
( |3.12| ). Suppose that we have a solution of the modified system ( |3.22| )-( [3.2S| ) with finite 
normahzation integral, ( ^.21 ). A direct calculation shows that the transformed functions 



a{r) = X-^a{X-\) , /3(r) = A"^ /3(A-V) (3.34) 

^{r) = A~^ ((/.(A-V) - v^oo) , w{r) = w{X~^r) (3.35) 



solve our original reciprocal limit system ( p.6[ ) , , ( |3.9D , ( |3.11 ) , and ( |3.12| ) with bound- 



ary conditions ( |3.13| )- (|3.15| ), if one replaces the energy E and coupling constant q by 

E = X-'^{E + ip^) and q = X^ q . (3.36) 



We point out that only the rescaled solutions ( |3.34| ),( ^35 ) and rescaled parameters ( |3.36| ) 



have a physical meaning. Therefore, we will in what follows consider only the rescaled 
tilde solutions; for ease in notation, the tilde will be omitted. 

In the remainder of this section, we describe our numerical solutions of the reciprocal 
limit equations. Just as in the case for the ED and EYM equations |^, there are 
solutions for different rotation numbers of the spinors in the (a, /?)-plane and for the YM 
potential in the {w, ii;')-plane. For simplicity, we restricted attention to solutions with 
rotation number zero for the spinors (as for the ground state solutions in [^). For the YM 
potential, we consider only the cases where the {w, t(7')-curve either makes a half rotation 
joining the points (1,0) and (—1,0), or makes a full rotation, ending at its starting point 
(1,0). A typical example for a solution of each type is shown in Figures || and ^. Because of 
the similarity of the YM potential to the BM ground state and the BM first excited state, 
we refer to these two types in what follows as the ground states and the first excited states, 
respectively. Notice that the curves in the (tu, t(7')-plane are not plotted all the way to their 
rest points at (1,0) or (0, —1), respectively. The reason is that we plot only the numerical 
solution on the interval (0,ri]. One sees that the spinors have become practically zero 
for r = ri, and it is thus an admissible approximation to smoothly join the (w, ti;')-curve 
with a vacuum YM solution by using the boundary conditions ( 3.33| ). We first discuss 



the ground state solutions. In Figure H the main characteristics of the solutions are 
plotted versus the coupling constant q. As explained above, E has the interpretation of 
the (appropriately scaled) energy of the Dirac particle. Since E is negative, the Dirac 



particle has gained energy by forming the bound state. The parameter B, (3.20), gives 
the total binding energy, i.e. the amount of energy which is set free when the binding is 
broken up. Since B is negative, we can expect that solutions of the full EDYM system, 
which are close to the solutions of the reciprocal coupling equations, should be stable. 
Finally, ryj and are the characteristic length scales of the solutions; more precisely, ryj 
is the radius where w changes sign, and Tq, is the radius where a has its maximum, 

w{ruj) = and o/{ra) = . (3.37) 

The characteristic radii are interesting because they give information about the "size" of 
the solutions as functions of r; i.e. they tell whether the fields are spread out in space, 
or whether they are localized close to the origin. It is also worth considering both radii 
because ryj and can behave quite differently (cf. Figure |3|). 

The plots in Figure | have a turning point at g ~ 8.49. Similar to the situation de- 
scribed for the spiral in ||], this is a bifurcation point which comes about as a consequence 
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of our rescalings. One branch of solutions can be extended up to q ~ 11.6. For solutions 
close to this end point, the potential w leaves the interval [—1,1] as shown in Figure ^ 
Since and both go to zero in this limit, the spinors and YM field are both confined 
to a smaller and smaller neighborhood of the origin. At the same time, the energy of the 
Dirac particle and the binding energy become infinite. The other branch of solutions ends 
near q = 8.95. For solutions near this end point, the {w,w')-cuTve comes very close to the 
origin before running into the rest point at (—1, 0), see Figure ^ This makes the numerics 
rather delicate, and we therefore have not yet analyzed this regime in much detail. It is 
interesting that is bounded near this end point, whereas r^„ seems to become infinite. 
This shows that, while the Dirac particle stays in a bounded region of space, the YM field 
becomes more and more spread out. 

For the first excited state, the energy spectrum and characteristic radii are shown in 
Figure ^. Since in general w never equals zero for the first excited state, we define via 
the minimum of w, i.e. 



w 



'(r^) = and a'ir^) = . (3.38) 



In contrast to the ground state, the solutions can now be extended up to g = 0. In 
this regime, the YM potential stays close to u; = 1; see Figure 0. The solutions have a 
bifurcation point at g ~ 9.866. The branch coming out at the bifurcation point for larger 
values of E is difficult to study numerically because the {w, w')-cuive comes close to the 
origin, see Figure ^. 

It is interesting that for the ground state solutions in Figure ^, the parameter q stays 
bounded away from zero, whereas the plots for the first excited state in Figure |^ could 
be extended up to g = 0. Let us consider how this can be understood directly from the 
equations. The parameter q enters only the YM equation ( |3.9[ ). In the limit q ^ 0, this 
equation goes over to the vacuum YM equation, which has only the trivial solution w = 1. 
Hence if we assume that the spinors have a finite limit for q ^ 0, then w{r) must go 
uniformly in r to one. This shows that the solutions can be regular for g — > only if 
w satisfies the boundary condition linir^oow{r) = +1. In particular, our ground state 
solutions cannot be regular in this limit. We next consider the limit g — > for the first 
excited state in more detail. Since w converges uniformly in r to one, the reciprocal limit 
equations (^^), (^), and ( 3.11| ) go over to the Dirac-Newton equations 



a' = -a - 2(3 (3.39) 
r 

/3' = {E + <f)a - -[3 (3.40) 
r 

r"^ ALp{r) = -c? . (3.41) 

These equations are obtained by taking the nonrelativistic limit of the ED equations |3| , 

and according to the results obtained in that paper, the equations (3.39)-(3.41) together 



with the normalization integral ( |3.12| ) have a countable number of solutions, characterized 
by the rotation number of the spinors (called the ground state, the first excited state, etc.). 
We thus expect that the functions (a, /3, ^p) corresponding to solutions of the reciprocal 
limit equations should for g — > 0, go over to a solution of the Dirac-Newton equations. The 
behavior of the YM potential w can now be analyzed in more detail by taking the solution 
(a, /3, ip) of the Dirac-Newton equations as a given inhomogeneity in the YM equation 
(|3.9D and performing a perturbation calculation for small q. More precisely, the ansatz 
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w{r) = 1 + q u{r) to first order in q, leads to the linear equation 



2 // 

r u 



2u + r a(3 



which can be solved by integration. Fixing the integration constants with our boundary 
conditions ti(0) = u(oo) = and n'(0) = 0, we obtain the unique solution 



u{r) 



ds 



s 



/ t a{t) f3{t) dt / ^ a{t) (3{t) dt 



1 



(3.42) 



This consideration shows that for (7 — > 0, the rotation number of w is uniquely determined 
by the rotation number of the spinors. Furthermore, one sees that in the limit g — 0, the 



Dirac wave function is determined by the Dirac-Newton equations ( 3.39 )-( 3.41 ). Thus only 
the gravitational attraction is responsible for the formation of the bound state, whereas 
the YM field has no influence on the spinors. 



4 Solutions of the EDYM Equations 

In this section, we shall construct numerical solutions of the full EDYM equations and 
discuss their properties. Our method is to first find special solutions which are small 
perturbations of either the BM solutions [Q] or solutions to the reciprocal limit equations of 
the previous section. We then trace these solutions while gradually changing the coupling 
constants. This yields one-parameter families of solutions which can be extended even 
to regions in parameter space where the solutions are far from all of the known limiting 
cases. 

In order to simplify the connection between the EDYM equations and the reciprocal 
limit equations of Section ^, it is useful to introduce a parameter e > in such a way that 
the reciprocal limit equations are obtained when e ^ 0. To this end, we parametrize the 
EDYM equations in terms of the new variables (e, g, E) as follows, 

n = (eg) 2 , e = 

1 

m = — — , oj = 

Since the EDYM equations involve three dimensionless parameters (namely m^K, uj/m, 
and e^), introducing {e,q,E) is merely a transformation to new independent parameters, 
prescribing at the same time the gravitational constant (this means that we give up the 
freedom to rescale r by fixing our length scale). In the limit e — > 0, both q and E go over 
to the corresponding parameters of the reciprocal limit system (see ( |3.8D and ( |3.5| )). Also, 
it is easy to check that the limits (|3)> ( CT) , and (|1|) are satisfied if we let e — > and 




keep (q, E) fixed. The parameters e and q can be written in dimensionless form as 



2 

m K 22 



e 



2 



m K e . (4.1) 



Thus e describes the relative strength of gravity versus the YM interaction, whereas q 
is the product of the gravitational and YM coupling constants. Up to a scale factor, 
E = uj — m. Since uj is the relativistic energy and m the rest mass of the Dirac particle, E 
can, exactly as in the previous section, be interpreted as the energy of the Dirac particle. 
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Finally, we also describe the binding energy by a parameter which corresponds to our 
notation for the reciprocal limit system ( 3.20| ) and set 



B 



p — m 



eq 



For the construction of numerical solutions, we use a two-parameter shooting algorithm 
combined with a rescaling method. Since this technique is quite similar to that described 
for the reciprocal limit equations in the previous section, we shall merely outline our 
procedure. In the first step of the construction, we consider the EDYM equations ( |2.23| ), 
(|]2|) and (|2.26D - (|2.28D with the side conditions 



< A" := 
< r := lim T(r) < oo 

r^oo 

lim w{r) = ±1 , 



(a^ + (3^) dr < oo 
VA 



(4.2) 

(4.3) 
(4.4) 



together with the following expansions near r = 0, 



a(r) 
A{r) 



air + 0{r^) 
1 + 0{r^) 



(3ir) 
T{t) 



1 + 0(r2 



For fixed e and g, we thus have the two parameters ai and E to characterize a solution 
of this modified EDYM system near the origin r = 0. On the other hand, we must satisfy 
two conditions at infinity; namely, w must converge to ±1, ( [4.4D , and the spinors must go 
asymptotically to zero in order for the normalization integral to be finite (|4.2| ). Hence, we 
can apply a two-parameter shooting method as described in Section ^. In order to have 
optimal boundary conditions at the upper end point r = ri, we again match this with the 
solution of the autonomous vacuum YM equation (see ( 3.32| ) and ( 3.33| )). The shooting 
method was again implemented in Mathematica, using an accuracy of 32 digits. For each 
solution constructed in this way, we verify that (^^) is satisfied and that the ADM mass 



is finite ( 2.32 ). Once we have found a solution of the modified equations, we rescale the 
solution according to 



a{r) = ^\-^a{\-'^r) 
A{r) = A(A~V) 



/3(r) = V^A-2/3(A-V) 
f{r) = T^^T{\-^r) , 



and transform the parameters (m, w, k, e^) as follows, 



m 
k 



A m 



X'e' . 



(4.5) 
(4.6) 



(4.7) 
(4.8) 



A straightforward calculation shows that the so-rescaled solution satisfies the EDYM equa- 
tions (|2?2^ ), (|]2|) and (|]2|)-(|]2D together with the original side conditions i ^lEj ) and 
( |2.3CI| )-( ^.33D . The parameters {e,q,E) transform under the rescalings as 



E 



A"^ {E + {t- 1) muo) 



(4.9) 



In the limit e — > 0, these transformations coincide with the rescalings of the reciprocal 
limit equations (|3.34| )- (|3.36| ). However, we remark that for the ED equations a much 
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different rescaling is used. Namely, in order to get a better correspondence to the reciprocal 
limit equations, we here scale the gravitational constant k, whereas in |Q k is fixed to be 
1 throughout. Clearly, only the rescaled solutions have physical significance. Therefore, 
in what follows we will consider only the rescaled solutions and again omit the tilde. 

In a realistic physical situation, the gravitational coupling is very weak, whereas the 
YM coupling constant is of order one, i.e. m?K <C 1 and ~ 1. Hence, according to ( |4.1| ), 
we will here only investigate the parameter range e <C 1, and we are particularly interested 
in the situation for small q. In the limit e ^ 0, there are known solutions of our EDYM 
system, namely the BM solutions (which more precisely are solutions in the limits k ^ 
and c^/k — > 1), and the solutions of the reciprocal limit system constructed in Section ^ 
We take these special solutions as the starting point for the numerics. By varying the 
parameters e and q and tracing the solutions with our shooting and rescaling methods, 
we obtain a two-parameter family of solutions. In order to reduce the computational 
workload, we did not step systematically through the two-parameter space, but always 
kept one parameter fixed while varying the other parameter. Since £ remains unchanged 
under the rescaling (see ( |4.9| )), it is most convenient to construct one-parameter families 
of solutions for different, fixed values of e. 

We now describe the solutions we found. Exactly as for the reciprocal limit system in 
Section ^ we restricted attention to solutions with rotation number zero for the spinors 
and a rotation angle of vr or 2tt for the YM potential. We again refer to these types of 
solutions as the "ground state" and the "first excited state," respectively. For the ground 



state solutions, the energy spectrum and the characteristic radii are in Figures U and 10 



plotted for different values of the parameter e (the characteristic radii are again defined 
by ( |3.37|) ). The curves A for e = coincide with the plots for the reciprocal limit system in 
Figure |3|. For small values of the parameter e, there are solutions for the EDYM equations 
which are close to the solutions of the reciprocal limit equations (compare the curves A 
and B in Figure |9|). In this parameter regime, the EDYM solutions look typically as shown 
in Figure 11. The metric functions A and T are both close to one; thus the gravitational 
interaction is weak, in agreement with our considerations after (3.4). The spinors and 
the YM potential look very similar to the solution of the reciprocal limit equations in 
Figure |5[ We conclude that the reciprocal limit system of Section |^ indeed describes a 
significant limiting case of the EDYM equations. However, one also sees that even for 
small e, not all the solutions of the EDYM equations are close to the reciprocal limit 
solutions. More precisely, curve B leaves the vicinity of curve A at g « 10 (see Figure |To|). 
If one follows curve B after it branches off from curve A, the parameter q first increases 
up to a turning point, and then decreases to g = 0. If e gets large, the solution curves no 
longer come so close to the reciprocal limit solutions (see curves C and D). The maximum 
of q decreases (see curve C) and finally disappears (see curve D). Figure [l^ shows a typical 
solution for small q. We note that in this parameter region, the metric functions A and 
T are not near one; this explains why the reciprocal limit equations are no longer a good 
approximation. Indeed, the potentials w, A, and T now resemble a BM solution of the 
EYM equations [|l| , and the spinors look like the solution of the Dirac equation in the BM 
background. Hence q ^ corresponds to the limit of weakly coupled spinors; i.e. spinors 
in a fixed BM background. Notice that the characteristic radii go to zero and the energies 
go to infinity in the limit g — > (see Figure P). This can be understood from our rescalings. 
Namely, for the (unsealed) solutions of our modified EDYM system, the BM solutions are 
easily obtained by taking the limit ai — > (in which the spinors go uniformly in r to zero). 



In this limit, the normalization integral (|4.2|) tends to zero, and thus the rescalings (4.5)- 
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( [4. 91 ) lead to a singular behavior of the rescaled solutions for g — > 0. To summarize, there 
is a one-parameter family of solutions (obtained by continuously changing the coupling 
constants), connecting the BM solutions to our reciprocal limit solutions 

We remark that our plots of the curves B have a small gap at q ^ 8.7. The reason 
is that in this region the numerics become unstable, and could not be carried out with 
our methods. But we were able to construct two branches of solutions which approach 
the problematic region from both sides. We suspect that the instability of the numerics 
is merely an artifact of our rescaling method, but it might well be an indication for a 
possible bifurcation point in this region. For the other curves C and D, we analyzed only 
the branch of solutions which extends towards smaller values of q. 

For the first excited state, the energy spectrum and characteristic radii are plotted 
in Figures 13 and ^ (the characteristic radii are again defined by ( ^) ). The curves 
A for e = correspond to the solutions of the reciprocal limit equations in Figure |6[ In 
contrast to the situation for the ground state, the solutions for small e are all close to the 
reciprocal limit solutions (compare the curves A and B). Figure ^ shows a typical solution 
for large q; one sees that the spinors and YM potential look similar to those in Figure ^ 
The form of the energy spectrum and the characteristic radii gradually change when e is 
increased; for example, the cusp in the ((?,r^)-plot becomes smooth (see curve D). It is 
interesting that for (7 — > 0, the curves converge independent of e to a single limit point 
(see Figure |l^. This limit point was already described at the end of Section ^ as the case 
when the spinors form a bound state due to their gravitational attraction, and the spinors 
generate a YM field (see ( 3.3g| )-( |3.42| )). This picture is in agreement with our numerics, 
since the spinors and metric functions, for a solution near this limit point, look similar to 
the ED solutions Q in the Newtonian limit, and ~ 1 (see Figure 16). The fact that 
this limit point is independent of e follows, because as explained at the end of Section ^ 
for g — s- 0, the YM equation decouples from the ED equations. For clarity, we point out 
that it would not be correct to say that the gravitational interaction dominates the YM 



interaction in the limit q ^ 0. Namely, according to (4J), the ratio of the gravitational 



and YM coupling constants is kept fixed, and thus g — > corresponds to the limit where 
both coupling constants go to zero at the same rate. Nevertheless, the YM field has for 
g — > no influence on the energy spectrum and the characteristic radii. 

A main qualitative difference between the ground state and the first excited state 
is that for the first excited state, we could not continuously join the solutions of the 
reciprocal limit equations with a BM solution. In order to see how this comes about, we 
did numerical calculations starting with a Dirac particle in the BM background (similar to 
that shown in Figure 17) and gradually increased the coupling of the spinors to gravity and 
to the YM field. For these "deformations of the first excited BM state," the curves of the 
energy spectrum and the characteristic radii have spirals, whose size and shape drastically 
changes when e is increased, see Figures IS and |l^. In the parameter regime where the 
energy plots spiral around, the spinors have self-intersections similar as observed for the 



ED solutions |3[, see Figure 20 



We now discuss the stability of our solutions. The relevant parameter for the stability 
analysis is the total binding energy B. Namely, if B is negative and smaller than the total 
energies of all other states, then energy is needed to break up the binding or to make 
a transition to any of the other states, and therefore for physical reasons the solution 
must be stable. Clearly, this energy argument does not provide a rigorous stability proof, 
and it also cannot replace the numerical analysis of linear stability (like e.g. in Q or |^ 
Section 8]), but it gives a strong indication for stability and is therefore commonly used 
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(see e.g. [^). Let us first apply this energy argument to the ground state solutions of 
Figures ^ and |l^. One sees that the total energy becomes negative for large q. For the 



curves B and C, this region is plotted in more detail in Figure 21. For the solutions on 



branch b, the total binding energy is minimal, and thus this branch is stable. Applying 
Conley index methods with q as the bifurcation parameter (see |1[1|]), we obtain, as in 
, that the two other branches a and c are unstable. Indeed, the instability of branch c 
follows also from the continuity of the Conley index and the fact that in the limit g — > 0, 
this branch goes over to the ground state BM solution which is known to be unstable 
[|[. When e is increased (see curve D, Figures ^ and 10), only one branch of solutions 
remains, which comprises the BM solutions as a limiting case and is therefore unstable. 
More precisely, the one-parameter family has in this case no bifurcation points, and in 
the limit g ^ the solutions tend to an unstable BM solution. Thus using Conley index 
techniques, it follows that the entire one-parameter family is unstable. We conclude that 
for small e, there is a stable branch of ground state solutions for which q lies in a finite 
interval away from q = 0; all other ground state solutions are unstable. For the stability 



of the first excited state, we consider the plots of Figures 13 and 14. Since for g — > 0, 



the spinors and metric functions go over to the Newtonian limit of the ED solutions, we 
conclude from |Q that the branch of solutions starting at q = should be stable. This is 
in agreement with our above energy argument, because on this branch the total binding 
energy B is negative, and is smaller than the total binding energy of the second branch of 
solutions, which comes out of the bifurcation point located at the maximum of q. Again, 
Conley index theory yields that this second branch is unstable. For the deformations of 
the first excited BM state, the total binding energy is positive (see Figures 19 and pO|), 
and hence these solutions should be unstable. Indeed, for the branch of solutions which 
extends up to g = (i.e. before the first bifurcation point), this also follows from the 
continuity of the Conley index and the instability of the first excited BM solution. 
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Figure 1: Reciprocal coupling limit: The ground state for q = 8.49811, E = 0.377446. 




Figure 2: Reciprocal coupling limit: The first excited state for q = 6.96132, E = 0.227762. 



E 




Figure 3: Reciprocal coupling limit: The energy spectrum and characteristic radii for the 
ground state 
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Figure 4: Reciprocal coupling limit: The ground state for q = 11.5838, E = 512161. 





Figure 5: Reciprocal coupling limit: The ground state for q = 8.76701, E = 0.334974. 
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Figure 6: Reciprocal coupling limit: The energy spectrum and characteristic radii for the 
first excited state 



20 




Figure 7: Reciprocal coupling limit: The first excited state for q = 0.285092, E = 0.16431. 




Figure 8: Reciprocal coupling limit: The first excited state for q = 9.75946, E = 0.405186. 




Figure 9: The energy spectrum and characteristic radii for the EDYM ground state and 
£ = (A), e = 0.0003933 (B), e = 0.005322 (C), and e = 0.02067 (D). 
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Figure 10: Details of the energy spectrum and characteristic radii for the EDYM ground 
state and e = (A), e = 0.0003933 (B), e = 0.005322 (C), and e = 0.02067 (D). 
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Figure 11: The EDYM ground state for q = 8.8373, E = 0.3367, and e = 0.0003933. 
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Figure 12: The EDYM ground state for q = 0.005881, E = 17867.5, and e = 0.005322. 
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Figure 13: The energy spectrum and cliaracteristic radii for tlie first excited EDYM state 
and £ = (A), £ = 0.0003933 (B), £ = 0.005322 (C), and e = 0.02067 (D). 
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Figure 14: Details of the energy spectrum and characteristic radii for the first excited 
EDYM state and £ = (A), £ = 0.0003933 (B), £ = 0.005322 (C), and £ = 0.02067 (D). 
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Figure 15: The first excited EDYM state for q = 9.548, E = 0.4153, and e = 0.0003933. 
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Figure 16: The first excited EDYM state for q = 0.09013, E = 0.1634, and e = 0.02067. 
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Figure 17: The deformation of the first excited BM state for q = 0.00047737, E = 367616, 
and £ = 0.02067. 
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Figure 18: The energy spectrum and characteristic radii for the deformation of the first 
excited BM state for e = 0.005322 (C) and e = 0.02067 (D). 




Figure 19: The energy spectrum and characteristic radii for the deformation of the first 
excited BM state for e = 0.0003933 (B),e = 0.005322 (C), and e = 0.02067 (D). 
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Figure 20: The deformation of the first excited BM state for q = 0.0005528, E = 729921, 
and £ = 0.0003933. 
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